#ifndef F8FF92F2_3E80_44E5_BC84_234F8A6413F1
#define F8FF92F2_3E80_44E5_BC84_234F8A6413F1

#include "functions.hpp"
#include "types.hpp"
#include <cmath>
#include <functional>
#include <math.h>

struct CalculationResult {
  double result;
  double error;
  int n;
};

class IntegralCalculator {
private:
  Func *func;
  IntegralMethod method_type;
  std::function<double(std::function<double(double)> const &, double, double,
                       int, double)>
      method;
  double a;
  double b;
  double epsilon;

  static double left_rectangles(std::function<double(double)> const &func,
                                double a, [[maybe_unused]] double _b, int n,
                                double h) {
    double total = 0.0;
    for (int i = 0; i < n; ++i) {
      auto v = func(a + i * h);
      if (isinf(v) || isnan(v)) {
        continue;
      }
      total += v;
    }
    return h * total;
  }

  static double right_rectangles(std::function<double(double)> const &func,
                                 double a, [[maybe_unused]] double _b, int n,
                                 double h) {
    double total = 0.0;
    for (int i = 1; i <= n; ++i) {
      auto v = func(a + i * h);
      if (isinf(v) || isnan(v)) {
        continue;
      }
      total += v;
    }
    return h * total;
  }

  static double middle_rectangles(std::function<double(double)> const &func,
                                  double a, [[maybe_unused]] double _b, int n,
                                  double h) {
    double total = 0.0;
    for (int i = 0; i < n; ++i) {
      auto v = func(a + (i + 0.5) * h);
      if (isinf(v) || isnan(v)) {
        continue;
      }
      total += v;
    }
    return h * total;
  }

  static double trapezoid(std::function<double(double)> const &func, double a,
                          double b, int n, double h) {
    double fa;
    double fb;
    bool fa_exists = true;
    bool fb_exists = true;

    // Check if f(a) and f(b) exist
    fa = func(a);
    if (isinf(fa) || isnan(fa)) {
      fa_exists = false;
      a += h; // Shift interval start
    }

    fb = func(b);
    if (isinf(fb) || isnan(fb)) {
      fb_exists = false;
      b -= h; // Shift interval end
    }

    double total = (fa_exists && fb_exists) ? (fa + fb) / 2 : 0;

    for (int i = 1; i < n; ++i) {
      double x = a + i * h;
      if (x >= b)
        break; // Stop if outside the interval

      double v = func(x);
      if (isinf(v) || isnan(v)) {
        continue;
      }

      total += v;
    }

    return h * total;
  }

  static double simpson(std::function<double(double)> const &func, double a,
                        double b, int n, double h) {
    double fa;
    double fb;
    bool fa_exists = true;
    bool fb_exists = true;

    // Check if f(a) and f(b) exist
    fa = func(a);
    if (isinf(fa) || isnan(fa)) {
      fa_exists = false;
      a += h; // Shift interval start
    }

    fb = func(b);
    if (isinf(fb) || isnan(fb)) {
      fb_exists = false;
      b -= h; // Shift interval end
    }

    double total = (fa_exists ? fa : 0) + (fb_exists ? fb : 0);

    for (int i = 1; i < n; ++i) {
      double x = a + i * h;
      if (x >= b)
        break; // Stop if outside the interval

      double coef = (i % 2 != 0) ? 4.0 : 2.0;
      double v = func(x);
      if (isinf(v) || isnan(v)) {
        continue;
      }

      total += coef * v;
    }

    return h / 3.0 * total;
  }

public:
  explicit IntegralCalculator(FunctionType f_type, IntegralMethod m_type,
                              double a, double b, double e)
      : method_type(m_type), a(a), b(b), epsilon(e) {
    switch (f_type) {
    case FunctionType::polynomial:
      func = Func::create_polynomial();
      break;
    case FunctionType::linear:
      func = Func::create_linear();
      break;
    case FunctionType::sinus:
      func = Func::create_sinus();
      break;
    case FunctionType::scaled_logistic:
      func = Func::create_scaled_logistic();
      break;
    case FunctionType::hyperbola:
      func = Func::create_hyperbola();
      break;
    case FunctionType::sqrt_hyperbola:
      func = Func::create_sqrt_hyperbola();
      break;
    default:
      func = Func::create_linear();
      std::cerr << "Invalid function type";
      exit(1);
    }

    switch (method_type) {
    case IntegralMethod::Trapezoid:
      method = IntegralCalculator::trapezoid;
      break;
    case IntegralMethod::Simpson:
      method = simpson;
      break;
    case IntegralMethod::LeftRectangle:
      method = left_rectangles;
      break;
    case IntegralMethod::RightRectangle:
      method = right_rectangles;
      break;
    case IntegralMethod::MiddleRectangle:
      method = middle_rectangles;
      break;
    }
  };

  CalculationResult calculate(int n) const {
    double result;
    double err;

    if (auto [fa, fb] = func->get_domain(); a < fa || b > fb) {
      std::cerr << "Invalid domain. Function is outside the domain\n";
      return {.result = NAN, .error = NAN, .n = n};
    }

    while (true) {
      double h = (b - a) / n;
      double coefficient;
      if (method_type == IntegralMethod::Simpson) {
        coefficient = 15;
      } else {
        coefficient = 3;
      }

      result = method(func->get_function(), a, b, n, h);
      err =
          std::fabs(result - method(func->get_function(), a, b, 2 * n, h / 2)) /
          coefficient;
      if (err < epsilon) {
        break;
      }
      n *= 2;
    }
    return {.result = result, .error = err, .n = n};
  }
};

#endif /* F8FF92F2_3E80_44E5_BC84_234F8A6413F1 */
